%% 函数或者脚本说明
%{
- --------------------------------------------------
* 文件名:ECM
* 函数名:None
* 功 能:通过等效电路方法计算吸收率
* 变量说明:
* 注意事项:
- --------------------------------------------------
%}
clc; clear; close all;
Frequency_range = linspace(.4, .8, 1000);

%% 数据读取区域
HTL_epsilon_real = readCST("HTL_epsilon'" + '.txt');
HTL_epsilon_real = HTL_epsilon_real(2:end - 1, :);
HTL_epsilon_imag = readCST("HTL_epsilon''" + '.txt');
HTL_epsilon_imag = HTL_epsilon_imag(2:end - 1, :);
HTL_epsilon_real_reshape = interp1(HTL_epsilon_real(:, 1), HTL_epsilon_real(:, 2), Frequency_range);
HTL_epsilon_imag_reshape = interp1(HTL_epsilon_imag(:, 1), HTL_epsilon_imag(:, 2), Frequency_range);

HTL_epsilon_real_reshape = HTL_epsilon_real_reshape + HTL_epsilon_imag_reshape * 1j;
clear HTL_epsilon_real;
% 读取CST计算的吸收率
Abp1 = readCST(".\type2\" + "Abp1" + '_type2.txt'); Abp1 = interp1(Abp1(:, 1), Abp1(:, 2), Frequency_range);
Abp2 = readCST(".\type2\" + "Abp2" + '_type2.txt'); Abp2 = interp1(Abp2(:, 1), Abp2(:, 2), Frequency_range);
Abp3 = readCST(".\type2\" + "Abp3" + '_type2.txt'); Abp3 = interp1(Abp3(:, 1), Abp3(:, 2), Frequency_range);
Abp4 = readCST(".\type2\" + "Abp4" + '_type2.txt'); Abp4 = interp1(Abp4(:, 1), Abp4(:, 2), Frequency_range);
Abp_broad = readCST(".\type2\" + "Abp_broad" + '_type2.txt'); Abp_broad = interp1(Abp_broad(:, 1), Abp_broad(:, 2), Frequency_range);

%% 定义一些式子
% 通用式子
c = 3e8; epsilon_0 = 8.854187817e-12; mu_relative = 1;
h1 = 30e-6; h2 = 20e-6; h3 = 50e-6; h4 = 35e-6; %介质层的高度
omega = Frequency_range * 2 * pi * 1e12;
Z_0 = sqrt(4 * pi * 1e-7 / epsilon_0); % 真空阻抗377
% 特用式子
epsilon_relative = HTL_epsilon_real_reshape;
beta = omega / c .* sqrt(mu_relative .* epsilon_relative); % 波数,只和材料相关
Z_c = beta ./ (omega .* epsilon_relative .* epsilon_0); % 材料的特征阻抗 sqrt(mu/epsilon)

%% 计算区域
Z_d = @(h) 1j * Z_c .* tan(beta .* h); % 定义不同高度介质的阻抗函数

% 求解Z_FSS
syms Absorption_CST_i Z_fss_i
syms Z_d_i
Z_total_i = 1 ./ (1 ./ Z_d_i + 1 ./ Z_fss_i);
Gamma = ((Z_total_i - Z_0) ./ (Z_total_i + Z_0)); % 反射系数
%{
% Absorption_ECM_i = 1 - ((Z_total_i - Z_0) ./ (Z_total_i + Z_0)).^2;
% Z_fss_answer_i = solve(Absorption_ECM_i - Absorption_CST_i, Z_fss_i); clear Absorption_ECM_i Z_total_i
Z_fss_answer_i = solve(Gamma - abs(sqrt(1 - Absorption_CST_i)), Z_fss_i); clear Absorption_ECM_i Z_total_i
% 符号计算得到的两个解
% - 1 / (Absorption_CST / (2 * Z_0 * (1 - Absorption_CST)^(1/2) - 2 * Z_0 + Absorption_CST * Z_0) + 1 / Z_d)
% 1 / (Absorption_CST / (2 * Z_0 + 2 * Z_0 * (1 - Absorption_CST)^(1/2) - Absorption_CST * Z_0) - 1 / Z_d)
% Gamma 符号计算
% - 1 / ((abs(Absorption_CST_i - 1)^(1/2) - 1) / (Z_0 + Z_0 * abs(Absorption_CST_i - 1)^(1/2)) + 1 / Z_d_i)
%}

Impedance = @(Absorption_CST_i, Z_d_i) (- 1 ./ ((abs(Absorption_CST_i - 1).^(1/2) - 1) ./ (Z_0 + Z_0 .* abs(Absorption_CST_i - 1).^(1/2)) + 1 ./ Z_d_i));

% 对于频率点1
% Absorption_CST_i = Abp1; Z_d_i = Z_d(h1);
% ans1 =- 1 ./ (Absorption_CST_i ./ (2 * Z_0 .* (1 - Absorption_CST_i).^(1/2) - 2 * Z_0 + Absorption_CST_i * Z_0) + 1 ./ Z_d_i)';

% ans2 = 1 ./ (Absorption_CST_i ./ (2 * Z_0 + 2 * Z_0 .* (1 - Absorption_CST_i).^(1/2) - Absorption_CST_i * Z_0) - 1 ./ Z_d_i)';
% impedance1 =- 1 ./ ((abs(Absorption_CST_i - 1).^(1/2) - 1) ./ (Z_0 + Z_0 .* abs(Absorption_CST_i - 1).^(1/2)) + 1 ./ Z_d_i);
impedance1 = Impedance(Abp1, Z_d(h1));
% 对于频率点2
impedance2 = Impedance(Abp2, Z_d(h2));
% 对于频率点3
impedance3 = Impedance(Abp3, Z_d(h3));
% 对于频率点4
impedance4 = Impedance(Abp4, Z_d(h4));

%% 测试吸收率
%{
Z_total_i = 1 ./ (1 ./ Z_d_i + 1 ./ impedance1);
plot(Frequency_range, real(impedance1), '-.', Frequency_range, imag(impedance1), ':');
Abs_ECM = 1 - abs((Z_total_i - Z_0) ./ (Z_total_i + Z_0)).^2;
figure;
plot(Frequency_range, Abs_ECM,'o', Frequency_range,Abp1,'*')
legend('ECM','Abp1');
figure; plot(Frequency_range, real(impedance1));
legend('R');
%}
%% 绘制阻抗
%{
% k1 = {'figure name', 'Impedance';
%     'title name', ' ';
%     'x_label', 'Frequency(THz)';
%     'y_label', 'Impedance(\Omega)';
%     };
%     k2 = {'x', 'y', 'legend', 'color', 'LineStyle';
%         Frequency_range, real(impedance1), 'real', 'b', '-';
%         Frequency_range, imag(impedance1), 'imag', [0.5, 0.5, 0], '-.';
%         };
YW_Plot({'figure name', 'Impedance';
    'title name', ' ';
    'x_label', 'Frequency(THz)';
    'y_label', 'Impedance Real(\Omega)';
    }, ...
    {'x', 'y', 'legend', 'color', 'LineStyle';
    Frequency_range, real(impedance1), 'real', 'b', '-';

    }, ...
    {
    'y_label2', 'Impedance Imag(\Omega)';
    'second y', 1;
    'second legend', 0;
    }, ...
    {'x', 'y', 'legend', 'color', 'LineStyle';
    Frequency_range, imag(impedance1), 'imag', 'R', '-.';
    });
%}

%% 并联在一起 第一种并联方法 XX不能用
%{
imp_parall = 1 ./ (1 ./ impedance1 + 1 ./ impedance2 + 1 ./ impedance3 + 1 ./ impedance4 + 1 ./ Z_d(h1));
% impedance_all = 1 ./ (1 ./ Z_d_i + 1 ./ Z_fss_i);
%}
%% 并联在一起 第二种并联方法
%{ open
imp_single1 = 1 ./ (1 ./ Z_d(h1) + 1 ./ impedance1);
imp_single2 = 1 ./ (1 ./ Z_d(h2) + 1 ./ impedance2);
imp_single3 = 1 ./ (1 ./ Z_d(h3) + 1 ./ impedance3);
imp_single4 = 1 ./ (1 ./ Z_d(h4) + 1 ./ impedance4);
imp_parall = 1 ./ (1 ./ imp_single1 + 1 ./ imp_single2 + 1 ./ imp_single3 + 1 ./ imp_single4);
%}
Absorption_parall = 1 - abs((imp_parall - Z_0) ./ (imp_parall + Z_0)).^2;
plot(Frequency_range, Abp_broad, '-.g', Frequency_range, Absorption_parall, 'r:');

% Out_port = [Frequency_range',Absorption_parall'];
% save Absorption_parall.txt -ascii Out_port
